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ABSTRACT 


A slender two-dimensional wedge wing at a high Mach number 
is subjected to an impulsive change in the angle of attack. The time 
history of the flow, till the final steady state is reached, is analysed 
and the unsteady flow solutions determining the temporal variations 
in the lift, drag, and rate of heat transfer are numerically obtained. 
The used parameters defining the problem are the free stream 
Mach number M., Reynolds number Re. * based on the free stream condi- 
tions and the length of the wedge L*, the semi-wedge angle 0 and the 
initial and final angles of attack oF and Ce. In the problem considered 


M., and Re, [ec 1 and the angles involved O° Ais Oe << 1. However, 


j 
their relative magnitudes are such that the combination Ney (eenge).”- 
remains finite and the product M,|9| < 0(1), where 6 is the maximum 
deflection of the wedge surfaces. 

Two types of thermal conditions on the wedge surface are 
used: insulated wedge surface and the wedge surface maintained at a 
constant temperature. The fluid is assumed to be a perfect gas with 


constant specific heat at constant pressure che constant ratio of the 


specific heats y, and constant Prandtl number Pr. 
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CHAPTER I 
STATEMENT OF THE PROBLEM AND ASSOCIATED LITERATURE 


1.1 Introduction 

Wings with the cross section of narrow wedge are commonly 
used as lifting surfaces in hypersonic vehicles. The behaviour of 
Such wings in steady flow conditions have been studied in great detail 
Over the past two decades. However, the question of maneuverability 
and control of hypersonic vehicles often necessitates the study of 
the flow in various unsteady situations. The present work deals with 
one such problem. A slender two dimensional wedge wing moving at a 
high Mach number is subjected to an impulsive change in the angle of 
attack. The time history of the flow, till the final steady state is 
reached, is analysed and the unsteady flow solutions are used to 
determine the temporal variations in the lift, drag, heat transfer 
rate and other characteristics of the wing. 

The disturbances due to the presence of the wedge are confined 
within a curved shock beyond which the free stream conditions prevail. 
For a slender wedge with a sharp leading edge the shock is attached to 
the leading edge. So the flows on the two sides of the wedge become 
independent of eacn other, which can be considered separately. Thus 
the problem under investigation reduces to the pone of finding the 
unsteady hypersonic flow over an inclined flat plate after an instan- 


taneous increase or decrease in the angle of inclination. 
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The hypersonic viscous flow over a slender body is characterized 
by the presence of a high temperature, low-density boundary layer adjacent 
to the body, which is considerably thicker than the corresponding boundary 
layers in subsonic or moderately supersonic flows. The displacement 
effects due to this thick boundary layer significantly changes the 
external inviscid flow and increases the pressure on the surface of the 
body. The flow in the boundary layer is, in turn, governed by the 
pressure (its streamwise gradient, in particular) and other external 
flow conditions. Thus the boundary layer and surrounding inviscid flow 
are mutually interdependent. This phenomenon which has been variously 
referred to in the literature as the hypersonic viscous interaction, 
viscid-inviscid interaction or shock-boundary layer interaction, is the 
basic feature that distinguishes the viscous flow at large Mach numbers 
from those at smaller Mach numbers. A review of literature on the 
interaction problem in two-dimensional hypersonic viscous flows over 


Slender bodies is given in the following Section. 


1.2 Review of Literature 

Shen [1,2] first considered the viscous effects in the steady 
hypersonic flow over slender wedges. He considered the entire flow 
region between the wedge surface and the leading edge shock as the 
boundary layer and established the validity of the boundary layer 
equations in this region provided the ratio of the boundary layer 
thickness to the axial distance is much less than unity. 

Detailed study of the hypersonic viscous flow by Lees and 


Probstein [3,4] helped to clarify some of the basic concepts. They 
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_ showed that the model postulated by Shen in which the boundary layer was 
assumed to occupy the entire region between the body surface and the 
shock, was inconsistant with the continuity condition. Figure 1.1 shows 
the flow regions on an inclined flat plate according to the model in- 
troduced by Lees and Probstein. The region between the shock and the 
body surface is divided into two distinct regions, the external layer 
and the boundary layer. In the external layer the effects of viscosity 
are negligible, so that the flow in this region is governed by the 
inviscid Euler equations. In the other region the flow is governed by 
the boundary layer equations (Shen [1], Lees [5]). Thus the hypersonic 


viscous interaction problem, mentioned earlier, has the following three 


aspects : 
SHOCK WAVE 
oy 
EXTERNAL 
LAYER 
—_ 
note SF 
_ BOUNDARY 


= LAYER 
tg Vda 
a 4 PLATE 


we Lor —S 
—>P> I / 


Figure 1.1 Flow Regions Over an Inclined Flat 
Plate in Hypersonic Flow 
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(i) The flow in the external layer can be considered as the in- 
viscid hypersonic flow over an “effective body" given by the 
boundary layer displacement thickness added to the original 
body. 

(ii) The distribution of pressure, tangential velocity and total 
enthalpy at the edge of the "effective body" determines the 
flow in the boundary layer. 

(iii) The displacement thickness used to obtain the external 
inviscid flow should be consistent with the displacement 
thickness obtained from the solution of the boundary layer 


equations. 


Various methods used to tackle the interaction problem for 
steady flow are reviewed in detail by Hayes and Probstein [6], Dorrance [7], 
Stewartson [8], Moore [9]. The basic characteristics of these approaches 
are outlined below. 

The flow in the external layer, which is equivalent to the 
inviscid hypersonic flow over a sharp-edged body of arbitrary shape can 
be obtained by the method of characteristics, the shock expansion theory 
and the tangent wedge approximation. The last method is the easiest of 
the three, gives results of sufficient accuracy and has been employed 
extensively in the hypersonic interaction problems. In the tangent 
wedge approximation the pressure at any point on a-slender body of 
thickness Oe is approximated by the pressure across an oblique 


shock that produces the local flow deflection dy%/dx* For free stream 
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Mach number M_ >> 1 and for a slender body with dy&/dx* <<eqyerfethe 
product K(x*) = M,, dy4/dx* remains finite, then the tangent wedge 


approximation gives the following results: 
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Linnel [10] obtained these results in a somewhat different form by 
assuming isentropic flow in the shock layer. Goldsworthy [11] later 
pointed out that the variation in entropy introduces errors of o(M.*) 


only. For small and large values of K equation (1.1) can be expanded 


as: 
p*(x*) a 2 
Sore adteiky ee ceo (1.4) 
and 
Balt e ail ; ; 
e- . 4 Ke ie at Lk an O(K 4) (1.5) 
PS Y (y+1) 


Equations (1.4) and (1.5) are often used in place of the complete 


equation (1.1) for the sake of simplification. 
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For the hypersonic interaction problem K in the equation (1.1), 


(1.4) or (1.5) has to be replaced by: 


d i * * 
Ub 4 d5*(x*) 5 (1.6) 


where yf (x*) defines the body surface and 6*(x*) is the boundary layer 
displacement thickness. Since the displacement thickness 6* in (1.6) 
is not known a priori, the use of tangent wedge approximation still 
retains the basic iterative nature of the interaction problem. How- 
ever, a great simplification is offered because the inviscid flow in 
the external layer need not be calculated explicitly. A disadvantage 
of using the tangent wedge approximation is that the vorticity in the 
external layer due to the curvature of the shock, which may affect the 
structure of the boundary layer, can not be accounted for. However, 
Hayes and Probstein [6] have pointed out that for a sharp-edged slender 
body the effects of the vorticity interaction can be neglected beyond 
a region very close to the leading edge. 

Lees and Probstein [3,4] showed that the extent of the inter- 
action between the boundary layer and the inviscid external flow is 


characterized by the hypersonic interaction parameter x , defined by 
X= MS (ee) (1.7) 
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peux 
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and C is the coefficient of the linear viscosity-temperature relation. 
Two extreme cases corresponding to y < 0(1) and y >> 1 were called the 
weak and strong interaction flows. The flow on a semi-infinite flat 
plate in the two extreme cases of weak and strong interaction has re- 
ceived considerable attention in the past. For the weak interaction 
case [3,4,6] corresponding to a region far from the leading edge, the 
appropriate expansion for the tangent-wedge approximation is given by 
equation (1.4). It was assumed that the pressure and other inviscid 
quantities varied slowly enough so that an approximate closed-form 
expression for 6* obtained from the zero-pressure gradient boundary 
layer solutions could be substituted in the expansion (1.4), resulting 
in a pressure distribution of the form p*/p* = 1 + Ay. A first order 
correction to the skin friction was then calculated as a perturbation 
to the zero-pressure gradient solution and it was shown that the heat 
transfer was unaltered to first order. 

In the strong interaction case, corresponding to a region 
near the leading edge, equation (1.5) is the appropriate expansion of 
the tangent wedge approximation. Using only the first term of equation 
(1.5) an order of magnitude analysis shows that in the strong inter- 
action region p* a y and 6*/x* a xe Seieaipeg Beef the 
boundary layer admits similarity solution. Solutions for this zero- 
order strong interaction theory were presented by Shen [2], Li and 
Nagamatsu [12], Stewartson [8,12] and Dorrance [7] for an insulated 
flat plate with Pr = 1. Li and Nagamatsu [14] later considered the 
case of a constant temperature plate. The constants of proportionality 


in the equations for p* and 6*, obtained by these authors were in general 
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agreement, the differences being mainly due to different methods of 
solving the boundary layer equations. In the works of Stewartson 
(and also Bush [15]) however, the tangent-wedge approximation was not 


“Ife and &* o, x#/4 the 


used. Instead, it was shown that for p* a x* 
flows in both the boundary layer and the external layer admitted similar 
solutions. A solution based on the method of inner and outer expansion 
was presented in which it was possible to take into account the effects 
of vorticity in the external flow. 

Lees and Probstein later generalized the strong and weak 
interaction theories for flow over a flat plate at non-zero angles 
of incidence [6]. The zero-order solutions are not affected and the 
effects of the non-zero angle of incidence appear through terms of 
higher order. In Section 3.2 the expansion scheme for the strong 
interaction case, outlined in Reference 6, is developed and solutions 
upto third order are presented. 

In the strong and weak interaction theories mentioned above, 
the pressure variations were such that the boundary layer equations 
admitted similarity solutions. The partial differential equations 
then could be transformed into ordinary differential Seas in 
suitable transformed variables. Such simplifications are not possible 
for moderately high values of xy. A theory valid for all values of x 
was presented by Nagakura and Naruse [16]. They considered the two- 
dimensional flow on an insulated slender body and rect the Karman- 
Pohlhausen integral method to solve the boundary layer flow for an 
arbitrary pressure distribution. However, in the final expression 


the pressure at any point on the body was given implicitly through a 
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9 
complicated nonlinear equation. This method is more suitable for the 
inverse problem of determining the body shape for a given pressure 
distribution than for the direct problem of determining the pressure 
distribution on a body of given shape. 

An approximate method, originally suggested by Lees [17] and 
later modified by Moore [18] is often used in interaction problems for 
general values of x. In this method, called the "local flat plate 
similarity" method the pressure gradient term in the momentum equation 
is neglected altogether, thus reducing the problem to the zero-pressure 
gradient Blassius problem. The basis of this assumption is the fact 
that the pressure gradient term in the momentum equation has a coefficient 
of the form y-1 which may become very small for the high temperature 
gas in the boundary layer. The solutions for slender wedges presented 
by Cheng et al [19], Dewey [20] and Mirels and Lewellen [21] belong 
to this category. Other approximate theories include the works of 
Bertram and Blackstock [22] and White [23] where a displacement thick- 
ness distribution of the form 6*/x* = axp*/(Mp*) was assumed. The 
coefficient a was found to remain approximately constant for the 
entire range of x, depending only on the thermal conditions on the 
body. The pressure distribution p*(x*) was then obtained from this 
relation and the tangent wedge equation (1.1). In spite of their 
approximate nature these theories often give quite useful and practical 
results. . 

Mann and Bradley [24] presented a different approach to the 
interaction problem which eliminated many of the approximations made 


in the references mentioned earlier. The flow over a flat plate at 
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zero incidence was considered. The non-similar boundary layer flow was 
obtained by numerically solving the partial differential equations. The 
external flow was solved by the method of characteristics. The two 
solutions were iterated till they converged within the required degree 
of accuracy. 

Two other solutions of the interaction problem by Kurzrock 
and Mates [25] and Butler [26] should be mentioned because of their novel 
approach. In both works the steady flow over a flat plate at zero 
incidence was obtained as the asymptotic limit of the unsteady "impulsive- 
Start-from-rest" problem. The complete unsteady Navier Stokes equations 
were solved numerically at discrete time steps over a region that ex- 
tended beyond the shock surface. Thus the interaction between the 
external layer and the boundary layer did not appear explicitly. The 
intermediate timesteps were considered essentially as iteration steps 
for the final steady solution and the authors did not attribute any 
significance to the unsteady solution. These methods require too much 
computer time to be of practical use in the general case. However, the 
results are useful to confirm the validity of the boundary layer assump- 
tions. 

The foregoing discussion was on the interaction problem in 
steady hypersonic flows. In unsteady flows the interaction problem is 
further complicated because the lateral velocity of the unknown displace- 
ment surface also affects the external flow. Lighthill [27] presented 
an extension of the tangent wedge approximation to obtain the pressure 
on a slender moving body of shape yor" st") placed in an inviscid hyper- 


sonic stream. This extension was based on the "piston analogy" of 
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Hayes [28], who observed the equivalence between the steady two dimensional 
Flow over an arbitrary body and the unsteady one-dimensional flow in 
front of a moving piston. The pressure in front of a piston moving 
with a velocity w* in a stationary compressible medium is given by 
2 iz 


* * * . 
Perey eee ge a = 4 (1.9) 


fey) 


where the subscript ~ refers to conditions in the undisturbed gas. 
Equations 1.1 and 1.9 are analogous if one identifies w* with the normal 
velocity component ua dy%/dx* induced because of the slope of the body. 
For a moving body an additional term oya/ot* is introduced in the normal 
velocity component. Thus the instantaneous pressure pa(x*, t*) on the 


body is given by equation ].1 with 
K(x®, te) = WED. tt (1.10) 


where 
oy* ove 


i Sea NTS Cr 11) 


Miles [29] has discussed the above expression and several 
other variations arising from the expansions (1.4) and (1.5), for the 
unsteady pressure in inviscid hypersonic flow. The first two terms 
of the expansion (1.4) lead to the acoustic approximation formula, 
while the piston theory of Lighthill [27] correspond to the first four 
terms of (1.4). 


For the unsteady interaction problem ve in (1.11) should be 
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replaced by (yx + A*), where A*(x*,t*) is the displacement thickness 
for the unsteady boundary layers, introduced by Moore and Ostrach [30]. 
Examples of unsteady interaction problem include the work 
of Rodkiewicz and Reshotko [31] who considered the unsteady weak inter- 
action flow on a semi-infinite flat plate at zero-incidence after the 
free stream Mach number was impulsively increased by a small amount 
(of 0(1%)). The pressure on the plate was assumed constant for calcu- 
lating the unsteady boundary layer and a small unsteady perturbation 
in the pressure was calculated from the boundary layer solutions. This 
approach is justified because of the small change involved in the 
problem. Rodkiewicz and co-workers [32-37] later considered several 
variations of this problem. 
Unsteady interaction problems involving lateral motion of the 
body are solely confined to cases of harmonic oscillations. King [38] 
and Orlik-Riickemann [39] considered the interaction on an oscillating 
slender wedge. These are examples of quasi-steady analysis, that is, 
the flow at any instant was assumed to be the steady flow corresponding 
to the conditions prevailing at that instant. Interaction problems 
involving impulsive lateral motion of the body are not available in 


literature. 


1.3. Description of the Present Work 

In the present work the tangent wedge approximation and its 
extension to the unsteady case, discussed in the last section, will be 
used to describe the interaction between the external flow and the 


boundary layer. The parameters defining the problem are the free stream 
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Mach number M. Reynolds number Re. L* based on the free stream con- 
ditions and the length of the wedge L*, the semi-wedge angle of and 
the initial and final angles of attack oF and Oe. In the problems 


av oF Lehi 


considered M, and Re. ioe 1 and the angles involved O2 Os 


However, their relative magnitudes are such that the combination 
M7 (Rey, x)!" remains finite and the product M |6| < 0(1), where 6 
is the maximum deflection of the wedge surfaces. The thermal condi- 
tions on the wedge surface also affect the problem. Two types of 
thermal condition are used in this work - insulated wedge surface and 
the wedge surface maintained at a constant temperature. The fluid is 
assumed to be a perfect gas with constant ch > y and Pr. However, the 
assumptions Pr = 1 and y > 1, used in some of the earlier references, 
will not be made here. 

In Chapter II, the governing equations and their transformations 
are presented. Chapter III is devoted to the problem of steady flow 
on an inclined flat plate. A four term series solution for the strong 
interaction region is described in Section 3.2. A more elaborate finite 
difference solution is presented in Section 3.3 and the results of the 
steady flow case are discussed in Section 3.4. In Chapter IV the un- 
steady flow on an inclined plate following an impulsive change in the 
angle of inclination is described. Initial conditions in time are de- 
rived from the steady solutions of Chapter III. The finite difference 
method of solution is described in detail. In Chapter V is indicated 
the utilization of the solutions on the two sides of the wedge yielding 


the unsteady characteristics of the two dimensional wedge wing. 
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CHAPTER II 
THE GOVERNING EQUATIONS 


2.1 Compressible Boundary Layer Equations for Two-Dimensional 
Unsteady Flow 
The governing equations for the flow in the boundary layer 
are given by [8]: 


Continuity equation: 


Q 


* 
nema tagee (CRUDE aoe (0%) © 0 (2.1) 
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Momentum equations: 
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Energy equation: 


3H oH* OHF ,*_scoptoys bBo oH* 
ox(Sa tut ae Smet Opa tae) 
] ‘) ; OUz 
“i (] - Dy) oy* (u*u 3y*) C2on) 


where H* = ole + ur2/2 is the total enthalpy. These equations, to- 
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gether with the equation of state: 


p* = px o* T* (2.5) 


and the viscosity-temperature relation: 


Mies lal) (2.6) 


define the complete set of equations governing the flow in the boundary 


layer. If we write 
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equations (2.1 - 2.6) take the non dimensional form: 
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eigen tapal: 26. 7 dlive, el Prel). 3 OU 
Re, L* 0 Br dy (gy) aaa ay). (2.10) 
peo) (Cane 
and u = u(T) (2.12) 
where g = (2-13) 
U*p*L* 
and Re, yx i (2018) 


is the Reynolds number based on the length of the wedge and the free 


stream conditions. 


The viscosity-temperature relation (2.12) is given by the 


Sutherland formula: 
te 
p= P/E So (2.15) 
T+ a 


where S* is a constant which has a value 110°K for air. The formula 
(2.15) is difficult to work with and a linear viscosity-temperature 


relation of the form 


Mesa as (2.16) 
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suggested by Chapman and Rubesin [40] has been extensively used for 
compressible boundary layer problems. The constant C is taken as 

1 /!(14S*/ TAY CT, +S4/T*) where Th is the non dimensional wall temper- 
ature. Thus the viscosity given by the approximation( 2.16) correspond 
exactly to the Sutherland value at the wall, where the effects of 


viscosity are most pronounced. 


2.2 Transformations of the Governing Equations 
The governing equations are now transformed to bring them to 
a form suitable for numerical solution. The first step involves the 


introduction of the stream function, wy, defined by 
ee (2.17) 


and the application of Dorodnitsyn-Howarth transformation which used 
in conjunction with the linear viscosity-temperature relation (2.16) 
eliminates o from the transformed equations. In this transformation 


the independent variables x, y, t are changed to X, y, t where 


(2.18a) 


x | 
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y= | o dy (2.18b) 
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For any function o(x,y,t) one obtains the following relations: 


a6 _ 9 ay 9 
so = 204 Sy 39 (2.19a) 
OX oy 
ae 9 = (2.19b) 
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pee (2.20) 
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Applying the relations (2.19) - (2.21) on equations (2.8) and (2.10) 
and making use of (2.9), (2.11) and (2.16) we obtain: 
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The next step of transformation involves expanding the y 
dimension so that the derivatives with respect to y are brought to 
normalized forms. In the usual boundary layer problems where the 
pressure is a known function, it is possible to eliminate the term 
Cp/Re., Lx: appearing in the equations (2.22) and(2.23), by suitably in- 
corporating it in the transformed variables. This leads to equations 
with coefficients of 0(1) which are convenient to handle numerically. 
But in the present problem the pressure is an unknown function to be 
determined in accordance with the interaction equation. In order to 
obtain normalized coefficients in the transformed equations we intro- 


duce the function p(x,t) given by: 


p(x,t) = Pbotd (2.24) 
Pox 
where Bes the zero-order strong interaction solution for the steady 
state pressure and Po is a constant to be determined in the next Chapter. 
The function p is of 0(1). 
The independent variables are now changed from x, y, t to 


x. We t where 
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Noting that X = nae rae where 


(2.26) 


is the interaction parameter evaluated at the end of the wedge, we can 


write 
n=Ayx lf (2.27) 
where 
Ree x 
A= 5 (yl (2.28) 
C Py x 
0 “L* 
is a constant. For any function 9(x,y,t) we have 
{6 = 2 0 oe (2.29a) 
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and use the relations (2.29), then the equations (2.22) and (2.23) 


transform to 
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The coordinates x and t have not been distorted in the trans- 
formations. This facilitates the retrieval of physically significant 


results from the solutions of the transformed equations. Thus we can 


write (2.31) and (2.32) in the form 
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It should be remembered that in the context of equations (2.33) and 


(2.34) 
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whereas, in the context of equations (2.7) - (2.11) 
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However, for variables such as p, which are independent of y(and hence 
independent of y and n) the partial derivative with respect to x 


and t are same in both senses. 


2.3 The Interaction Equation 
For the unsteady flow over a flat plate inclined at a small 


angle 6, to the main flow, the shape of the "effective body" is given 


by 


Veet) Be AE(XAGT®) (2.35) 


a 
7 


ic 
bas (£€.8) enoiveups to txe% it Ss 
. ‘4 ~ 
i a 
7 ee i> 0s, a 
<@ —s ¥ ¢ 
leet 


lee 
(PFS) = (Ws) enot suns Yo axainog art n 
: we a 
1 ee Ke 
a 


» 
“~ 

re 

are 

i? ih Due ¥ 
| : ue 
Me, 
bith 
“aan *. | W 


sone bas)y to sehinelait ame dotdw a 98 Rue a 5 
x ot tosqesay Atiw oviteyinsb aor: om 


23 
From equation (1.11) the normal velocity component at the edge of this 


effective body is 


dA* 
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(2.36) 
Here A*(x*,t*) is the displacement thickness for the unsteady boundary 
layer, introduced by Moore and Ostrach [30]. In terms of the present 


notations the differential equation for A* becomes 
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Using the continuity equation for the external flow 


aoe (p*u*) = 0 (2.38) 
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If we apply the transformations of the last section on the two integrals 


appearing in (2.39) it is found that 
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Here terms of o(M~*) have been neglected from the integrand of the 
last integral. Writing 6* for both the integrals in (2.39) and using 


(2.38) once again, we obtain 
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A further simplification is obtained if ua 1s replaced=by us. “This 
is consistent with the observation of Stewartson [8] and Goldsworthy 


[11] that 
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u, = <= = 1 + :0(M) (2.42) 


when M, is large but the product of M, and the streamline deflection 
in the external flow remains finite. Therefore from (1.1), (1.11) and 


(2.35) we have 
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where 
” 26 
K(x,t) = M,(o, + $2 + 22) (2.44) 
and 
wai co | 
3/4 2 
6* = * 
S(x.t) =te= Eh th Bh er (2.45) 
YPo a P fe) 


The interaction problem can now be stated as follows: the 
functions f(x,n, t) and H(x,n,t) are to be solved from equation (2.33) 
and (2.34) for a function p(x,t) which is related to f(x,n,t) and 


H(x,n,t) through equations (2.43-45). 


2.4 Boundary and Initial Conditions 

The conditions on the wedge surface and at the outer edge 
of the boundary layer provide the boundary conditions in n. On the 
wedge surface both components of velocity are zero and the wedge is 


either insulated or held at a constant temperature. Thus, 


at y" = yF 
u* = 0 (2.46a) 
v* = 0 (2.46b) 
= 0 (for insulated wedge surface) (2.46c) 
or 
T* = TE = constant (for constant temperature wedge surface) (2.46d) 
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In terms of the transformed variables we have 


at n = 0 
nt On, (2.47a) 
of 
ite Um, (2.47b) 
oH 
on = 0 (2.47c) 
or H = H. = constant (2.47d) 


é At the edge of the boundary layer the tangential 


where Hy = ck TH/> ux 
velocity and the total enthalpy in the boundary layer should match with 


the corresponding quantities in the external layer. Thus 


aby ee 
u* = us (2. 48a) 
H* = He (2.48b) 


In terms of transformed variables we have 


at n> 
ot 1 (2.49a) 
Tonk (2.49b) 


where terms of order um? have been neglected. 
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In boundary layer problems it is customary to satisfy the 
conditions (2.49a,b) at a finite value n = Mee The value of We is 
chosen sufficiently large, so that, in the resulting solution 
a" F(n,)/ an" | Hin S)i2\end |a“H(ng)/an‘| » k > 1. can be made less 
than a predefined tolerance. From preliminary experiments it was 
found that, for the problems considered in the present work Nae 8 
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was sufficient for a tolerance of 10 ~. 


The initial conditions in x and t are discussed in Chapters 


til and IV. 
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CHAPTER III 
THE STEADY FLOW ON AN INCLINED FLAT PLATE 


3.1 Introduction 

In this Chapter the steady hypersonic flow on a flat plate 
inclined at a small angle 6 to the free stream will be considered. 
The solution of the steady problem will provide the initial conditions 
for the unsteady problem at t = O and also the final solution to which 
the unsteady flow should converge as t > ~, 

Specializing the governing equations to the steady case, we 


have from equations (2.33) and (2.34) 


3 2 2 2 
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For the steady problem p = p(x) is given by (2.43) with 


K(x) = M,(0, + $8) (3.3) 
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and 6 = 6(x) given by equation (2.45). 


3.2 Series Expansion Solution of the Steady Problem* 

A series expansion scheme for solving the steady strong 
interaction flow on an inclined flat plate, originally suggested by 
Lees and Probstein [6], is considered in this Section. Using appro- 
priate expansions for p, 6, f and H in terms of the interaction para- 
meter x and K. (= MO.) » the partial differential equations (3.1) and 
(3.2) are reduced to a sequence of ordinary differential equations. 
The unknown constants appearing in the expansions of p, 6, etc. are 
then determined from the solutions of these ordinary differential 
equations in accordance with the tangent wedge equation. 

The zero-order term of this scheme correspond to the strong 
interaction similarity solution considered in references [6,7,8,12, 
13,14] where the assumption Pr = 1 was made to simplify the energy 
equation. Solutions for higher order terms are not generally avail- 
able. In this section solutions for terms up to third order are 
carried out for arbitrary values of Pr. 


We start with the expansion 
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of the tangent wedge equation (1.1). If we define the constant 


“contents of this Section were presented in the Third Canadian Congress 


of Applied Mechanics, Calgary, 1971 and published in AIAA Journal, 
Vol. 9, No. 3, pp. 535-537, (1971). 
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K. 7 M84, (3.4) 


then from (3.3) we have 


Se Max (375) 


For the strong interaction region (x >> 1) an order of magni- 


tude analysis [6] shows that 


Ma 1/2 
do 
or pew (3.6) 


Writing K(x) = K, + constant qe and substituting in equation (1.5), 


we have for y >> 1 
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xX xX 


where p,(i=0,1,2,3) are constants to be determined. Consistent with 


(3.7) we have for 6 


where 6. (i=0,1,2,3) are constants. On substituting (3.7) and (3.8) 


in (3.5) and (1.5) one obtains the following relations between the two 
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sets of constants: 


Po = az Yiytl) 6 (3.9a) 
8 
P) = 3 (5) + 3) (3.9b) 
p = 10 6 fe 32 3yt] ] (3 9c) 
oi 3 EY, 
y(yt1)" 5 
se eats aA (3.9d 
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In view of the expansions (3.7) and (3.8), Lees and Prob- 


stein [9] suggested the following expansions for f and H: 
FUXsny = fy i9)) as ee, 


and 


Hi (n)K.@ @H,(7) +9 (aq) K 
H(x,n) = Ha (n i" areas Se Bie 


On inserting the expansions (3.7), (3.10) and (3.11) in the 
equations (3.1) and (3.2) and equating the terms of the same order on 


the two sides we get the following sets of coupled ordinary differential 


equations for f.(n), H.(n) (i=0,1,2,3): 
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Here the prime represents differentiation with respect to 
n- In equations (3.13) through (3.15) the unknown constants p,(i=1,2,3) 
appear in the inhomogeneous parts. To eliminate these unknowns we 


introduce a further change of variables: 
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= z\,2 g cate 
fz = (f3-P3f)/P]” > H3 = (Ha-P3H5)/p, (3.18) 


From equations (3.13) through (3.15) we obtain the following 
differential equations for f.(n), H, (n) (i=1,2,3): 


LY) + oh = FL (n) (3.19a) 
(2H) + 103) (Fy = FL2)(n) (3.19b) 


The linear operators iJ) (i,j=1,2,3) are given by 
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2 : 
F| )(y) = PrHof for, ‘t= <2 (3.25) 


+-2(1-Pr)L(Fh") + FFE] for i=3 (3.26) 


with Cy = 1.5 and Cy = Ze 


From (2.47) and (2.49) the boundary conditions can be 


written as 
f (0) = 0 (3.27a) 
F5(0) = 0 (3.27) 
Hg (0) = 0 (for insulated surface) (Sez7c) 
or Hy (0) = He (for constant temperature surface) (3.27d) 
fi(n,) = 1 f2o2ie} 
Ho(n@) = | ST) 


and for f,, H, (i=13253) 


F.(0) = 0 (3.28a) 
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F, (0) = 0 (3.28b) 

H, (0) = 0 (for insulated surface) (3.28c) 

or H. (0) = 0 (for constant temperature surface) "62280 
Fy (n,Q) = 0 (3.28e) 

Hs (ng) = 0 (3.28f) 


In order to determine the unknown constants Ps» 6; (i=0,1,2,3) 
we insert the expansions (3.7), (3.10) and (3.11) in (2.45) to obtain 


the following expression for 6: 
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Using (3.9) and (3.31) we now obtain the following ex- 


pressions for the constants P: 
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Numerical method of solving the ordinary differential 


equations for fy, Ho; f., H.(i=1,2,3) is described in the following 


subsection. From these solutions the integrals Ips I), In, I, are 
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calculated. The constants pss 6, (i=0,1,2,3) are then obtained from the 


equations (3.31) and (3.32). 


3.2a Method of Solution 
The equations for Footy are nonlinear, while those for f, 
H. (i=1,2,3) are linear. The problems are of the boundary value type, 
that is, three conditions are specified at n = C0 and two conditions 
at y= Nes 
Let us first consider the zero-order problem 


vue it 1c — 
inane fof’ + B(Hy-f ) = 0 (3.12a) 


vt i Eve Leite 
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Earlier workers often used the assumption Pr = 1, which eliminates 
the last term in (3.12b). For Pr = 1 and the insulated surface case, 
(0) = 0 (also for constant surface-temperature case with Hee 1) 
equation (3.12b) has the constant solution Hg (n) = ]. The zero order 


problem in such cases reduces to the classical Falkner-Skan problem: 
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(Op. 04; fj(n,) = 1 (3.33b) 


In their work on similar solutions of compressible boundary 
layers, Cohen and Reshotko [41] have studied extensively problems of the 
type (3.12,27) for arbitrary values of 8. It was shown that for 
positive values of 8, the conditions (3.27) are sufficient to give a 
unique solution. In the present case 8 = (y-1)/y is positive. 

The nonlinear boundary value problem given by equations 
(3.12,27) can only be solved iteratively. In the method employed by 
Cohen and Reshotko, the differential equations were transformed to an 
integral system with fo as the independent variable. At each iteration 
Stage a better approximation to the solutions was obtained by numerically 
integrating two functions that involved the current approximate solu- 
tions. In the present work we use a more direct approach, called the 
"shooting metnod", where the solution to the boundary value problem is 
obtained by solving a sequence of initial value problems. 

We first consider the insulated surface case. The equations 


(3.12) can be solved as an initial value problem if all the elements 


of the row vector 


are known. We first assume suitable approximations.o, and hy for the 


unknown conditions fy‘ (0) and Hy (0) . With the initial condition 
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the equations (3.12) are solved up to n = ne by the predictor- 


corrector method described in Appendix A. From the solution we note 
eng )wee] (3.36a) 
and B, = Ho (na) - | (3.36b) 


which represent the errors in satisfying the conditions at n = Nes 


The equations are solved two more times with the initial conditions: 


Cy = bO. 0: O, +e, » 0] (3374) 


[0, 0,0, , A, +, 0] (3.37b) 


and C3 
where « is a small quantity. Let (Ay Bo) and (Az .B3) be errors corres- 
ponding to (3.36) in the two cases. Interpolating from the above re- 


sults we obtain more accurate approximations for the missing initial 


conditions: 


e(B, A, - BA, ) 


% = 9 * R(By-By) + AplB,-B,) + AglB,-B,) ee) 
c 1 A, B.,-B. # A, By B, 7 Az B, B., 


e(ByA, - ByAp) 


ho = Ay + R(BocBo) @ A(BocE.) * ALB.<B.) 
2 1 A, B.-B3) 3 A,(B3-B,) + A3(B,-B.) 


(3.38b) 
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The calculations are repeated with these values of the missing conditions. 
The iteration is stopped when |A,| and |B,| becomes less than a prede- 
fined tolerence. The initial approximations om and dy are found by 
preliminary trials. The perturbation term e is progressively reduced 
as the solution approaches convergence. 

For constant surface temperature case the missing initial 
conditions are Fo‘ (0) and Ho (0) - The above iteration scheme is used 


with 


Cy = [0, Q, O1s Hi» hyd (3.39a) 
Co = [0, 0, O7 76. Hy i] (3.39b) 
Cz = POR 0, Oy» Hy» , + € | (3.39c) 


The problems for f. H.(i=1,2,3) are linear and therefore, 
can be solved directly. For each i we solve the following three 


initial value problems using the predictor-corrector method described 


in Appendix A: 
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HW) (9) = 0 s HO (9) = 0 (3.40) 


62) (9) = 0  Hf4)(0) = 0 (3.41) 


H63) (9) = 1; H89)(0) = 0 (for insulated surface) 
(3) form of (84 ten 
or H‘*/(0) = 0 ; H‘°’(0) = 1. (for constant-temperature surface) (3.42) 
The solution can then be obtained from 


= e(1) i n(2) 4 pe (3) : (3.43a) 
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where the constants A, B are adjusted to satisfy the conditions 3 (ng) 


and Hs (n,) = 0. A and B are given by 
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A= - < 3.44 
Fy Hn) - #2) HH +) are 


yo i i 3.44b 
37 ay; ( ) 


The FORTRAN program for the series solution is described in 


Appendix A. 


3.3 Finite Difference Solution for the Steady Problem 
The series expansion solution discussed in Section 3.2 has 


a 4 This solution is therefore suitable for 


an error of order x 
sufficiently high values of x only. In the problems considered in 
the present work Xp* can be as low as 0(1), which requires a more 


accurate solution. In the current Section a finite difference method 


is presented which refines the steady state series expansion solution. 


The complete form of the tangent wedge equation is used. 


For the finite difference solution it is necessary to specify 


the initial values of f(x, on); H(x, »n) and p(x, ) on some initial line 
X= xX). In numerical solutions of the steady boundary layer problems 


[22 ,42-44] the usual practice is to specify these initial conditions 


at x = 0. Since all the x-derivatives are multiplied by x, the partial 


differential equations (3.1) and (3.2) reduce to ordinary differential 
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equations in n at x = 0. The solutions of these ordinary differential 
equations are then used as the initial values. This approach is 
questionable because it assumes the validity of the boundary layer 
equations at x = 0. In the present work an alternate approach in 
line with that used by Baxter and Flligge-Lotz [45] is employed. The 
values given by the series expansion solution on a line x = X, near 
the leading edge are used as the initial values. Xy is given a 
value of the order TOT which corresponds to x ~ 0(10°) , so that the 
series expansion solutions are of sufficient accuracy at x = Xa 
Lees [5] mentioned (o*/x*)® Passa as the limit of the applicability 


of the boundary layer equations. The present choice meets this criterion. 


The initial conditions are, therefore, 


PyF Cn) Ky PoFp(n)+L5F,(n n)+p, “Fa (n) 1k, 


f(x, on) = fo ah a ye (3.45a) 
X] x ; 
pyHy(n) Ky poHy(n)+Lp3H,(n)+p,°H3(n) Jk, 
H(x, sn) = Ho(n) + tit _Palaln)+LegRn)+0y Fn)" (3.45b) 
X] x] 
p,K Potp. K 
p(x,) = 1+ is + 43} (3.45c) 


where xy = x* po and the other quantities are as defined in the 


previous Section. 


The interaction between the boundary layer and the external 


layer is tackled by an iterative procedure. Each iteration cycle con- 


Sists of the following two parts: 
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(i) to solve equations (3.1,2), for any given p(x) distribution, 
in the region xt x ane. 13:0) <tin <n, Subject to the initial 
conditions (3.45a,b) and boundary conditions (2.47) and (2.49), 

(ii) to obtain a new p(x) distribution from the solutions f(x,n), 
H(x,n) of (i), making use of equations (2.43),(2.45), (3.3) 
and (s..45¢). 
The steps (i) and (ii) are repeated until p(x) converges 

within a specified tolerance. Details of the two parts are given in 


the sub-sections 3.3a and 3.3b. 


3.3a Finite Difference Solution of the Boundary Layer Equations 

From the series expansion solutions it can be observed that 
the x-derivatives of different flow variables f, H, p, etc. are large 
for small values of x and the x-dependence decreases in the downstream 
direction. For the same degree of accuracy a finite difference scheme 
should therefore employ smaller stepsize near the leading edge. In 
the present work we use the following arrangement of x-stations, in 


which the stepsizes increase in the positive x-direction in geometric 


progression: 
AX, = X (3.46a) 


AX; = x - Xs] = K AX sy Acari ree 5 (3.46b) 


Here KY is a constant greater than 1 and J is an integer such that 
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ag Te 
In equations (3.1) and (3.2) the x-derivatives appear in 
the form 4x 3¢/dx where ¢ is a dummy function. A finite difference 


formula for this term for the present arrangement of x-stations can 


be written as: 
PO ae (b Ciao. aces 258 3 
[4x ax Bs. (B. Ci)o. 4 Cabs o> FEZ Sd (3.47) 
where o, represents o(x,) and the coefficients Bs C. are given by: 


B. = PURER sadiee 


J 
2 a oH P31=34 ited (3.48) 
C, = 0 , j=2 
: vs me Wei hee (3.49) 


This corresponds to three-point backward difference formula for j=3,4, 

...,J0 and two-point backward difference formula for j=2*. The trunca- 
: Z Z 

tion error involved in equation (3.47) is of the order x Ax; (0 6/9x ); 


for j=2 and x Ax5(3°9/ 3x") , for j=3745c ace. 


t j j j i ; 1 to zero 
For j=2, $3_9,is undefined. But since Co has been set equa 
ae es present any problem. This arbitraryness has been intro- 
duced for the sake of uniformity at all x-stations, j=2,3,...,d. 


at | | : * a 
. yitiwora ae 

nt Yasaqs ere sat (8.2) bas 4h.8) enottauparal 
gonsistttb aataet A .nottonut wyimub s 2f > srariw solos bm 

neo 2hotdete-x to sneepnsyIS: trsearq old ot arist etd aa 

als | :28 nest 


ee 


ee ahah “i oy : 


:vd nevip 946 2 te 8 2tnatoittees oft bag (x4 zinsesnaor | 


| | i ae 
Sf . ANC 8 = (8 


2 “i 
| (f+, 48)8 ear 
‘. | L,.--—h, & = * : . 
(ae £) y ¢xA cm | ‘ihe 
td = 2? | 7 
S$ 
x ye - 
(eh.€) ee > ae xi T Ne an 


-sonust sfdT .*S= yo? slumyot sonsistttb brswiosd tntog-ows ay : 
Cases 6) ,Xdpx ebro ada Yo et (SALE) nota sups nt bevfovar 108 


, g 
be enahsGxt 10? af a 


| : ~~ ie 
.>,€=t 10t sfumrot sonevettib biswissd intoq-sers of 2bnoge: 102 id 


O1SS OF a ts2 ia esd eae 
“orsat need zeonyyersidis 2 ? Nir 
devsalsSHt sonotontoes 1152 ia a 


47 
Substituting (3.47) in (3.1,2) we obtain 


D-fi''.¢ f)'T(14B.)f. - (BtC.)¢. 
pf, ; [(1 B.)f, (B.+C fs) + C.F. 5] 


grey 
ef i Be tOr) fo. = CBC.) f! .f! = 
GLB +05) f) RBG.) et C5 f5_o] + QjH, = 0 (3.50) 
and 
co . 
veetoH! +B.)f. - na Oe me if2 
en) HL Bt (Bits) T iy i Cif 5-24 
- f'lB.H. - .+C.)H. «He 
fi[B.H. (B. C5 )Hs + CH; 1 
2 
— (Pr-l 
i meats), 1 guys 
cer Tamla ul re ie (3.51) 


In equations (3.50,51) the prime represents differentiation 


with respect to n and the function Q(x) is defined by 


Q(x) = g[1 - ox Py (3.52) 
p 
At any station x = Xie eee Hey: Fi_o and Hs _2 are all known functions 
of n and Pys Oy. Bs» C, are known constants. Thus equations (3.50,51) 


form a set of coupled nonlinear ordinary differential equations for 


ur and H, of the form: 
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OgH' tag fH'+),(n)H'ta7F'H+O4(n)Fitag((f!')+FF'''] = 0 (3.54) 


where the a's are known constants and the ¢'s are known functions of 


n. The boundary conditions are: 


(0) = 0 
f'(0) = 0 
H'(0) = 0 
ie (3.55) 
H(0) = Hy 
thingler tl 
Ch ie 


By solving the ordinary differential equation problem 
(3.53-55) at x = Xoo Xgo vee oXy One can generate a numerical solution 
of the partial differential equations (3.1,2) in the region xX SX < Xj; 
On SA This approach, often called the "difference-differential 
method" in the literature, was originally suggested by Hartree and 
Womersley [46] and has been used extensively by Mann and Bradley [24] 
and Smith and co-workers [42,43,45] in boundary layer problems. Be- 
cause of the use of backward difference formulas for the x-derivatives 


this implicit method is inherently stable [42]. 
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Equations (3.53-55) can be considered as a fifth order non- 
linear boundary value problem with three conditions specified at n = 0 
and the other two specified at n = Ne: Such a problem can only be 
solved iteratively. The method most commonly used is the "shooting 
method" described in Section 3.2 in connection with the problem for 
zero-order functions fos Ho: 

The shooting method, although adequate for certain problems, 
has a severe limitation. In this method it is implicitly assumed 
that if the conditions at one boundary are perturbed by a small amount, 
the changes in the solution of the initial value problem near the other 
boundary will be of similar magnitude. But the equations (3.53,54) 
do not have such a behaviour. The coefficients of some terms in these 
equations involve the quantities B; and C; which are multiples of the 
term Xi /AX;. The stepsizes Ax. are to be kept small to reduce the 
truncation errors of the finite difference scheme. This makes the 
term X5/AX 5, and therefore some of the coefficients of (3.53,54), 
quite large as the solution proceeds in the downstream direction. 
The presence of these large coefficients makes the solution extremely 
sensitive to small changes in the initial conditions [42]. While 
experimenting with the shooting method, the present author found situ- 
ations where a change of 0(107°) in f''(0) could cause a change of 
0(10°) in F'(n,)- The shooting method is obviously unsuitable in such 
cases. Cebeci and Keller [47] introduced a modification called "parallel 
shooting" which reduced this sensitivity in certain cases of the 


Falkner-Skan problem. However, for the problem under consideration, 
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this modification did not offer much advantage. 

An alternative approach based on the "quasilinearization 
method" was developed for solving the problem (3.53-55). The quasi- 
linearization method, which is essentially an extension of the Newton- 
Raphson method of solving algebraic or transcendental equations (Bellman 
and Kalaba [48]), provides a powerful tool for solving nonlinear boundary 
value problems. Radbill [49,50], Libby and Chen [51] and Jaffee and 
Thomas [52] applied this method to equations arising from the boundary 
layer problems. 

In the present case we formulate an iterative scheme, each 
cycle of which involves solving linearized forms of equations (3.53,54). 
These linear equations are solved by a finite difference method where 
the conditions at both boundaries could be used simultaneously. Bell- 
man and Kalaba [48] and Kenneth and McGill [53] gave proof of con- 
vergence for second order boundary value problems. Similar analysis 
for the problems considered here are not available in the literature. 
However, at each iteration cycle we have to calculate two quantities 
which indicate how closely the current solution satisfies the original 
nonlinear equations. Examination of these terms indicate a uniform 
rapid convergence and also provides a suitable means of stopping the 
iteration after the necessary accuracy has been achieved. 

At each station x we solve for a sequence of functions 
fl] (ny, Hen) , m= 1,2,... starting with a suitable initial assump- 
tion flO], | LOT) , where each set of functions fim), HEM on), 


m= 0,1,2,... satisfies the boundary conditions (3.55). At any station 


spsinevbs doum Y9tto Jon bib nofss: 

norsestrsont frasup” ens mo voasd dasingas svi toned 6 Hh 
-t2sup SAT Age -£2.£) metderq oxtsesgatabiny tah begoleveb.2sw | 
-nodwah ant to ier oneets: ns vlletinezes ef aotriw ,bodsom notsss 
rvs f T38) enoitsupe fetnabnssensit +o otsidepis prtvioe to bodiem 
visbnuod weontfinon paftyvioe 10% food [utvswoq 6 2z9bfvor”d (£88) « 
bra ssttsel bas [fe] ned? bus yddid .(Oeveh] i tdbsh — | 
yi sbnuod ari movt pifers 2norssups at bord am ets bari ggs {s2}, es 


neg, satsnoe oveSs1est we ots Momo? ew. a289 sne2ex9g BAS nh) = 


= 


(D2 £218) anortsups to 2ntot pasiiseatt onivioe 2oviovat AoFriw toe 
siolw borktom saderettib eiratt a Ye boylo2 s%S 2nottaups “son Bi 
-[fod- .yfedenedfumte beev odebives esfasbmiod nied vs anots ‘ANOS 
-io2 to Fo6yd Svép- [be] Titaom brs rigenniod baa [8h] sdeNsd t bas: 
et ey! ‘sis xeltmre...emeidorg avlsy yibbaved abo brssee ini 
autsisstl off nt slfdsftsvs gon 9v6 ‘oned: borabtenos or 
exist inoup ows Sts usta os eo ow id spt sevad? hey) bie: 
sntptvo sid eottattee notiutoe Sneqs. aad Nis2ofs wots 
motiny s sisatbat amass seodt to fotisntwexd .znottoups 9 
ont Srna Gee +0 208M" eldsd rue. 6 asbivorg oels brs sr 
sbaysitrion, al zed Yoptuods wiseesodn afd isthe Notts 
eng Se to sonsupaz 6 “git aviez ow Ee nottste dosg TA. ih 
-qiwees istting aldetive 6 dat ontsiste vo vetel edo Cn My 
(a) bm a) hart: Yo doer tiase ovartw | mye) oleh Mi it 
noftsse ye tA phen £) anoistbnoo.yssbmiod oft eitetons « ven 5 O=n 
va 


- 5 baat 5 ae 
: eres ; if 7 ine ; 
y ae 


> «tei one 


5] 
x; the solutions of the preceeding station X;_1 are used as the initial 


assumption. Thus we set 
fl) =F, 4(n) (3.56a) 
Helm) = He y(n) (3.56b) 
If we introduce the functions e'")(n) and g(")(n) such that 
el] (ny = fl] (ny - ght (,) (3.57a) 


gl" (1) = pln) - qm 1], ,) (3.57b) 


then at the m th iteration cycle the momentum equation (3.50) gives 


[m]. 


the following linearized equation for e 


c eld: t ‘+0, [m](,)_bmd: ‘405 ET(n) el 45,5 b"I(n) ebm ae rim], 


(oyae) 
with the boundary conditions 
alM(g) = 0 (3.59a) 
el™ (9) = 9 (3. 59b) 
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and 


el™]i(y ) = 0 (3.59c) 


The coefficients in equation (3.58) are defined below: 


iS =p. 


: (3.60a) 


oN(n) = (14B,)F""H(n) - (By4C,)F, y(n) + C5F; p(n) (3-606) 


ogl™(n) = - 2(B,40,) 401 (n)+(By4C,)F1_4(n) - CyFt_p(n) (3.60¢) 


0g8™U(n) = (rap, el een) (3.604) 


and 


Fund (ny = py fen) 
+ elm-1V6n) (148, ye (np (B40), 4) + Csf,_9(n)} 


: FLT Ton) (B40) Tm) (BAC 5) fn ee 


+Q, yb Vn) (3.60e) 
The function plm](,,) represents the departure from satisfying the 


momentum equation (3.50) at the m th iteration cycle. The problem 


given by equations (3.58-60) is solved by a matrix method described 
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in detail in Appendix B. The solution el™(,,) is then used to obtain 
the next iterate fm] (,,) = fim-1])(,) + elM (,) 

From the energy equation (3.51) we obtain the following 


(m) 


linearized equation for g* ’: 


ghd + gl (nygh™l® + ol] (nygltd = - el™ (ny (3.61) 


with the boundary conditions: 


gl™J(9) = 0 for insulated surface (3.26a) 
or gl™](9) = Q for constant-temperature surface (3.62b) 
and gin) = 0 (3.62¢) 


The coefficients appearing in equation (3.61) are defined below: 
=p. .6 
ees pale? 


o4t"(n) = (14B,) 4" (n) - (B,#C,)#;_,(n) + C;f5 p(n) (3.636) 


o,f (n) = - B ft" (n) (3.63c) 
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(n) + C.f. ,(n)} 


+ Hem Tiny (148, £0 (n )-(B,+C;)f, j'j-20" 


isle 


- mln) ca Hl I, )=(B.+C.)H 


Te _(n) Hors 54y- o(n)} 


+ op, PED ccglNiny)? + ellen el (n)3 (3.634) 


Here again, the function glm(,,) represents the departure from satis- 
fying the energy equation (3.51) at the m th iteration cycle. The 
problem defined by equations (3.61-63) are solved by a matrix method 
im) ,) 


details of which are described in Appendix B. The solution g 


Hen) = Hem") ¢) + gl (1). 


is then used to obtain the next iterate 


The iteration is stopped after the m th cycle if 
maxt|fel™D] |p fg4™y]) prem) iets < (3.64) 
where ¢ is a predefined tolerance and ||9$|| for any function 9(n) is 


[Io] | = max f{[(n)]} (3.65) 
O<n<n, 


Once the criterion (3.64) is satisfied the functions pln] (.,) 


and yb] (.,) are accepted as the solution at the station x = Xj Hess 


f.(n) = FE (n) : (3. 66a) 


H.(n) = Hem (1) (3.66b) 
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This completes the calculations at the j th station and the solution 
proceeds to the next station Xeaye 


For a tolerance e = (le. the number of iterations necessary 


was 2 in most cases and was never more than 4. 


3.3b Solution of the Interaction Equation 
From the solution of the boundary layer equations the pressure 
distribution can now be obtained using the equations (2.43), (2.45) and 


(3.3). We first introduce the normalized quantities 6(x) and I(x): 


B(x) = M, A x 3/4 (3.67) 
SoX,* 
- a 
By ‘} (H-f'")dn 
I(x)*s Si ay ee (3.68) 
0 
So and Ip are constants associated with the zero-order solution de- 


fined in the Section 3.2. In view of the relation (3.3la) between the 


constants 5p, Po and Ig, equation (2.45) becomes 


p(x) = ar (3.69) 
6 (X 


From (2.43) we recall 
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This relation can be inverted to give the following expression for K 


in terms of p: 


(P Pg X- 1) 


K(x) = — cane 
Ly" + (Bpox - 1) MOFg1/2 


Since x = Xie riggs equation (3.70) can be written as 


where 


53 4 Po XL* 


and p has been replaced by I/é from (3.69). 
Now, from (3.3) 


3.-70)) 


(30) 
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(3.73) 


Writing K, = M8, and using the relations (3.67) and (3.72) we obtain 
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Equation (3.74) gives the following differential equation for 6: 
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The method of obtaining the pressure distribution p(x) from 
the solution of the boundary layer equations can now be described as 
follows. 

From the solution of the boundary layer equations f(x,n), 
H(x,n) we calculate the I(x) distribution using equation (3.68). 
Equations (3.76) and (3.71) then defines a nonlinear first order 
ordinary differential equation for 6 with the initial condition: 


at xX = Xy 


I P(x) 


oS ctx) = (sav) 
where f(x, sn); H(x, sn) and p(x) are those given in (3.45). We then 
solve for 6(x) and obtain the pressure function p(x) from (3.69). 

We solve the differential equation for 6 by an iterative 
method based on quasilinearization. Let aed (x), ny ="O8 128 . eae beta 
sequence of iterates satisfying the initial condition (3.77). If we 


define a sequence of functions 
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n 
then = I(x) represents the correction to be added to the existing solu- 
: n-] 
tion gt diss From (3.71, 76,77) we obtain the following linearized 
problem for cee 


x Fob x) gtr) - - oly) (3.79) 


ina) = 0 (3.80) 
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[n-1] 
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4 96 
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a, | gl Pe et lS 
cy x2 x) git aes xl/2 
and 
1/4 
= Keex [n-1] 
Ci cane | dé. 3~ s(x), 
Dv’(x) = ix = tF6 - + (3.81b) 
X 4 Ga SA | 


The superscript [n-1] in the last terms of (3.8la,b) implies that the 
current iterate gin-1] is to be used in place of 6 . The function 
pind) represents the departure from satisfying the differential 
equation at the end of the (n-1)th iteration cycle. ~hn](,) is 
solved at the stations Xs j=2,3,...0 by means of a finite difference 
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method described in Appendix C. The next iterate 
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is then calculated. The iteration is stopped when 
max {{[cEMl), | oldiy3 < « 
Where € is a predefined tolerance and for any function (x) 


[lol] = max {]o;|} . (3.82) 
1sjsd 
Once the criterion (3.82) is satisfied, gin is accepted as 


p. I U 6 2 >) j 2S wa J Sees 


which gives the pressure distribution corresponding to the I distri- 


bution. 


3.4 Results of the Steady Flow Problem 

The results of the series solution method of Section 3.1 and 
3.2 are presented in Tables 3.7 and 3.2 for Pr = 1.0 and Pr = 0.72 for 
y = 1.4. Corresponding results available in the literature are mostly 
confined for the zero-order solution and for Pr = 1.0. Comparison 
with the available results are given in Table 3.3. 

The lack of good agreement with the results of reference 
[13] is due to the fact that in Stewartson's work the tangent wedge 


approximation was not used. References [14] and [16] made use of 
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TABLE 33:1 
RESULTS OF THE SERIES SOLUTION, °Pr°="1.0, 7 = 1.4 


Pe ee 
Constant Temperature Surface 


= * = = = = 
Hy 1.0 Hy 0.75 H. 0.5 Hy 0.25 Hy 0.0 


#*"(0) 0.76275 0.70985 0.65555 0.59966 0.54192 
f, ‘(0) -0.43902 -0.38851 -0.33604 -0.28122 -0.22352 
f,'(0) -0.50238 -0.43317 -0.36097 -0.28517 -0.20488 
f3' (0) 0.37567 0.32262 0.26704 0.20833 0.14559 
H9 (0) 0.0 0.12767 0.25152 0.37118 0.48618 
H! (0) 0.0 -0.049415 -0.095230 -0.13684 -0.17343 
H5(0) 0.0 -0.044431 -0.083600 -0.11648 -0.14164 
H3(0) 0.0 0.030704 0.057669 0.080050 0.096598 
I, 1.3091 1.0802 0.84660 0.60775 0.36291 
I, 0.60845 0.49842 0.38816 0.27797 0.16835 
I, 0.59464 0.48268 6.37203 0.26383 0.15905 
I, -0.15318 -0.12451 -0.094934  -0.067762 -0.040542 
89 0.73394 0.66669 0.59022 0.50008 0.38643 
6 -0.80117 -0.88426 -1.0011 -1.1824 -1.5227 
85 -0.82398 -1.0034 -1.2861 -1.7975 -3.0026 
63 0.22395 0.27017 0.34326 0.47659 0.79750 
Po 0.50904 0.42008 0.32920 0.23632 0.14112 
Py 1.4969 1.6418 1.8486 2.1792 2.8402 
P. 1.5098 1.8139 2.2948 3.1765 5.3451 
P3 1.3067 1.5745 1.9986 2.7758 4.6750 


ee. OO 


“For Pr = 1.0, the constant-temperature surface case with He ="1, 0715 
identical with the insulated surface case. 
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TABLE 3.2 
RESULTS OF THE SERIES SOLUTION, Pr = 0.72, y = 1.4 


SR ee ee eee 
Constant-temperature Surface 
Insulated ©§ ——______ 


* = = = = = 
Surface Hh 1.0 H. 0.75 Hy 0.5 H. O#250)H 0.0 


bret) t$00. 73347 0.77184 0.71433 0.65514 0.59403 0.53065 
f, ' (0) -0.40862 -0,44635 -0.39195 ..-0.33525. -0.27576  -0.21274 
f5'(0) -0.45843 -0.50958 -0.43573 -0.35844 -0.27690 -0.18994 
fu). 0.33730 0.38013 0.32374 0.26443 0.20143 0.13354 


Ho (0) 0.83244 -0.078114 0.037865 0.15023 0.25859 0.36245 
H, (0) 0.0082451 0.036826 -0.0096634 -0.052783 -0.091917 -0.12620 
H4(0) 0.012646 0.039139 -0.0040486 -0.042388 -0.074853 -0.099956 
H3(0) -0.024620 -0.044235 -0.010735 0.019043 0.044208 0.063401 
Ty e275 13903 Tet3t6 0.86725 0.59655 0.31850 
I, 0.56600 0.6418] Bee WAS 5373509 pirate )es} 0.14609 
I, 0.54984 0.62411 0.49885 037543 0.25400 0.13710 
I, -0.17254 -0.17970 -0.14593 -0.11225  -0.079041 -0.047019 
59 0.70780 0.75637 0.68237 0259737. 0.49544 0.36201 
64 -0.83070 -0.77930 -0.86640 -0.99244 -1.1985 -].6318 
55 -0.88748 -0.77852 -0.96141 -1.2609 -1.8406 -3.4377 


53 0.21841 0.20040 0.24380 0.31388 0.44639 0.78795 


0.47342 0.54063 0.44001 Wess 8) 0.23196 OU1Z2385 


Pg 
Py 1.5524 1.4475 1.5976 1.8175 2.1863 3.0146 
Po 1.6183 1.4126 1.7194 2.2220 3.2053 6.0359 


P3 1.3305 1.1918 1.4507 rosa 2.6830 4.8985 
"For the insulated surface case the numbers given in the rows 5-8 are 


Hy (0) » H,(0), H, (0) and H.(0) respectively. 
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TABLE 3.3 
COMPARISON WITH THE AVAILABLE RESULTS (Pr=1.0, y=1.4) 


pa M00 (0) = Oy 


) 
Sou 
aos (H.=1.0) H,=1.0 -O0 (H,=0.0) (H)=0.0) (H,=0.0) 


a SS ER alle. ee a NC ee hr 
Stewartson [13] 0.555 0.881 


Li et al [14] 0.514 0.766 0.149 0.539 0.408 
Nagakura et al [16] 0.510 1.49 0.764 
Present work 0.50904 1.5098 0.76275 0.14112 0.54192 0.48618 


i 


the tangent wedge approximation. However, the method of solving the 
boundary layer equations were different in the two works. Li and Nagamatsu 
used analog computers while Nagakura and Naruse used von Karman's integral 


method. In references [13] and [16] a two term expression of the form 
Earels x + constant 


was obtained for the pressure on a flat plate at zero incidence. The 
constant term corresponds to PoPo in the present notation. 

In Figures 3.1 to 3.8 the results of the series solution are 
compared with those of the finite difference solution. The finite 
difference solutions were carried out with the first station Xt 1./1024 
and Ky = 1.15. This corresponds to 36 x-stations between x = 0.0 and 


x = 1.0. The step sizes in the n-direction were taken as 1./16 which 
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63 
corresponds to 129 n-stations. 
Figures 3.1 and 3.2 show the velocity and total enthalpy 
profiles at four different locations x = 0.25, 0.5, 0.75 and 1.0 along 
a flat plate at an angle 2° to the free stream with M thecOF0, X* = 8.0, 
Pr = 0.72, y = 1.4 and Hi = 0.5. Since the error involved in the series 


solution are of order pelle 


the agreement improves in the upstream 
direction. The overall agreement is good because of the high value 

of Xp* . Figures 3.3 and 3.4 show the variation of the surface shear 
function foo. the surface enthalpy gradient H,;, the nondimensional 
pressure and displacement thickness p and 6 along the plate obtained 
by the series solution and the finite difference solution. The series 
solution predicts lower values for Fis Hh and p and higher values for 


6, the discrepancies being of the order of 10%, 5%, 15% and 10% at 


The agreement between the series solution and the finite 
difference solution decreases with the decreasing values of Xp * be- 
cause of the larger error terms in the series solution. The disagree- 
ment becomes significant for Xp* <4. In Figures 3.5 - 3.8 the re- 
sults from the two solutions are compared for a flat plate inclined 
at an angle of 2° to the free stream with M, = 10.0, Pr = 0.72, 

y = 1.4 and Hh = 0.5 for Xp * = 2.0. The disagreement is quite notice- 
able as can be seen from the velocity and total enthalpy profiles in 
Figures 3.5 and 3.6. The variations of foo Hes p and 6 along the 


plate are shown in Figures 3.7 and 3.8. The values predicted by the 
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series solution at x = 1.0 are about 70%, 20% and 25% less for foe 
Hy and p respectively while the value for 6 is about 50% higher. 

In the remaining figures of this section the results from 
the finite difference solutions are presented to show the effects of 
changing the parameters b? His X* on foo Hy » Prand 6 . 

Figures 3.9 through 3.12 shows the distribution of Fees 
6 and p along the plate for 6, ee? Oey 2 ode, O° OWT P Yee Os rc. 

y = 1.4, M = 20, H, = 0.5 and Xp * = 8. 

Figures 3.13 through 3.16 presents the effect of increasing 
Xp* on the distribution of fi', Hp s 6 and p along the plate for 
prea On) ey = beAS NS TOs here and 0, = alee 

Figure 3.17 through 3.20 shows the distribution of Foe Hh 
6 and p along the plate for Hp = 0.0, 0.5, 1.0 with Pr = 0.72, y = 1.4, 


M,, = 20, Xj» = 8.0 and o = 2°. 
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CHAPTER IV 
THE UNSTEADY PROBLEM 


4.1 The Initial Conditions 

In this chapter we consider the unsteady viscous flow on a 
plate, whose inclination to the free stream has been impulsively changed 
to OD Ff at t = 0 from the original inclination he at. t.< 0. ™ The 
steady state solution of the last Chapter, corresponding to Oe hae 
describes the flow on the plate prevailing at t < 0. At the instant 
of the impulsive change (t=0), the spatial distribution of the velocity 
and temperature field with respect to the flat plate is given by this 
steady state solution. However, the change in the inclination of the 
plate changes the shape of the effective body which in turn changes 
the pressure distribution. 

The initial conditions in x are specified at x = xy for 
t > 0. In accordance with the approach in the steady problem, we 
assume that the series solution corresponding to oe = Poet gives the 
solution at x = x, for t > 0. This is equivalent to assuming that the 
final steady state is reached instantaneously at x = xy This assump- 
tion is justified because the close proximity to the leading edge makes 
tne] TOW dbex = xy extremely insensitive to the angle of inclination 
and the initial and final steady state solutions there differ by 


less than 1%. Thus the initial conditions in x are given by 
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fOr t>.0 
F,(n)P)Ky ¢ — Pofy(n)+[paFy(n)+p,-F_(n) IK- 
feces) fy(n) + its Se EER TE) 5 
xX) Xy 
H, (n)p4K pH, (n)+[p oH. o(n)+p,°H 3(n) 1K 
H(x, nt) = HA(n) + lin Oa , Paelnl+Legfgin)+Py "FNS 6 4 9 
] 0 /2 
1 X 
p, K Pptp. K 
p(x, .t) = 1+ ety 2-3 bf (4.3) 
1 Xx] 
—- _- -1/2 2 
where mM hte and Ky of eit On of: 
The initial conditions in t are given by 
tor. XxX. xy 
f(xsn,0) = #6) (xyn) (4.4) 
H(xsn.0) = H' (xn) (4.5) 


where (0) and 469) are the steady state solutions corresponding to 
8, = OH i? obtained by the finite difference method of Section 3.3. 
The initial pressure distribution p(x,0) = 519) (x) is calculated by the 
method of Section 3.3b using Ky of in place of K,- The initial and 


boundary conditions are shown schematically in Figure 4.1. 


4.2 Finite Difference Method for the Unsteady Problem 

The unsteady equations are solved at discrete time-steps 
t., 1=1],2,3,... with stepsizes increasing progressively in the positive 
t-direction: 
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Figure 4.1 
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on ABEH: f = 0; "= Oe H= "Beton ="0 
on FCDG: f'= 1; H=1 
on ABCD: £ = f°; H=H°; p= p® 


on ADGH: f= f,; H= Hy; P= Py 


Arrangement of the Boundary Conditions (Schematic) 
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At, (4.6a) 


ty - ty] = Aty = Katy) » i-2,3,... (4. 6b) 
Here the first timestep At, is of 0(1079) and K, is a constant greater 
than 1. The arrangement of x-stations is the same as that given in 
(3.46). Figure 4.2 shows schematically the arrangements of x and t 


stations. 


Figure 4.2 Schematic Arrangement of the x and t Stations 
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We introduce the notation 


FY (sn) = f(xsnsty) (4.72) 
HO) (xn) = H(xansty) (4.7b) 
pil (x) = p(x,t,) ete. (4.7c) 


For any function ¢(t) the first derivative in the backward difference 


form can be written as: 


(2) ae uot?) - (u-tv, gl) is y,o0te) (4.8) 


where the quantities u., V; are given by 


K. +1 
] t 
WQS ere for i = | 
1 At. Ky 
1 2K, +1 Mm 
ares aes pop ae ct (4.9) 
pee 
Woe 0 for i = | 
2 
ale coe for i > 2 (4.10) 
«At, OKyt = ; 


The truncation error in equation (4.8) is of the order of at, (a@o/at?) 
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for i=1 and At.*(a%¢/at®) (4) for 4 > 2. 


Using (4.8) equations (2.33) and (2.34) can be written as 


ac 3 2 2 2 
Wiest 4 a ax(2t (2 4 y,) - Eat 


3n° dn ‘dxon Ui OX an 
(i-1) (i-2) : 
taxt(ujty) Se = vy SE xy tHE = 0 (4.1) 


, 2(Pr-1 mt int (2) FLARE (4.12) 


Here absence of any superscript in f and H indicates Fi) (yn) and 


ae g(x) and ri) (x) are given by 


6) 
Q(x) = aft - Fey By (4.13) 
p 
Tt) 
and abi) (x) = 4g a (SP) (4.14) 


The use of backward difference formulas for the t-derivatives 
results in an implicit scheme which is inherently stable. At each 
time step t = t. i1=1,2,3,... we approximate the functions ee 
gf) and pli) by extrapolating from the values at the previous time- 
steps. The partial differential equations (4.11,12) are then solved 


in the region xX; <x < Xj); OF le ss using the initial conditions 
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(4.1-3) and the boundary conditions (2.47,49). The method described 


in Section 3.3a is followed. From the solution we then obtain pe 
in accordance with the interaction equations (2.43-45). The solution 


then proceeds to the next time step t = t. The solution is continued 


Llp 
till ap/ot becomes less than 10° in magnitude and is then compared with 
the steady state solution corresponding to 6 = On Ff The method of 
solution for the time steps ty and ty is different from that for the 
general time step t., i > 3. They are described in the following sub- 


sections. 


Le) 


4.2a General Case, t = ty es 


At the time steps t. (i>3) we use the following extrapolation 
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where D. 
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z, = - 48 aT (4.18a) 
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2K, +2K,+1 
z= —+——y (4.18b) 


(2K, +1) Ky 


be 


- 7 aan 
bodiozeb boriiam SAT . (28, We.S) anoitttbnos yxsboued sit bis (eh 
(i> pistdo madd ow notsufoe ot mova Bewohfot 2t a 


he 


a 
a - 


aoftuboe ofT .(@A8).5) enotisups eo ont Ajtw sons 
bountsnos ef notsulo2 sat py -d =a J ee sai xen ont oF ebssoong 
dgiw bevsqnioo vedt 2f bas obustrpenm nf “oF oe 229 zomored so\ge Tf 
to bodsom ont y, 4° = § oF ntbiaqesyid9 nattufee stste- | aa 
ant vot s6Add mort snaystttb Sf gf bas ,J 2qste - ott Ot not tu 
-duz pniwolfot oft nf badirozeb ots yenT  .e, < i ey at qote omits Ts 


E< ta, = $ . 926) levensd 


notisfeqs1sxs putwof fot sit seu sw (E<t), -J eqase airs ott tA 
(i, as Hp | (Hy vot esfum 


tay (Dey ott tguy egylge > ae 


as) . to) Ypgee (oP tpn = 9p 


expe Pgeeay Wg seta 


(487.8) Sar giS) +3 é i 5 ms 
oat 
—— os mi 
‘ Ja 
: 7 7 Bae by 


92 
2 


ae. (2K, “+K, +1) (K, +1) 
al Wee SEES a a (4.18c) 
(2K,+1)K, 
1+K, 
he en 1g (4.18d) 
t 


The truncation errors in equations (4.15), (4.16) and (4.17) 
are of the order of at.?(a@p/at?) 1), at,2(a@qvat?) 44) and at f (329/802) 
respectively. 

With these values for pil), gti) and R‘) the equations (4.11) 
and (4.12) can now be solved in the region Aiesehes XGs On on: 
The method discussed in Section 3.3a is followed. The x-derivatives 
are replaced by backward difference formulas and the resulting ordinary 
differential equations are solved at the stations x ;(J=2,d) using the 
iterative method based on quasilinearization. The linearized problems 
to be solved at the m th iteration cycle are given by (3.58-63) with 
the following exceptions: 
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- ax Cu, Fl Hin)-(ujav, 4, 1-1) (n) + eta) (4.19b) 


on) =. B fl) - 4xsus +R, (i) (4.20a) 
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Here Lm (,.) and ylm(,) are the m th iterate for 1) (n) 
and HS1) (n). The extra terms in (4.19,20) over the corresponding ex- 
pressions (3.60c,e) and (3.63c,d) result from the t-derivatives in 
(2.33,34). These extra terms do not present any additional difficulty 
and the matrix method described in Appendix B is used to solve the 
linearized equations for elm] and gin], 

After solving IG) and 561 (n) at all x-stations 
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TS $14 7" j=2,9 (4.21) 
J Io 
The pressure distribution pi) can now be obtained from the equation 
(2.43-45) by the method described in Section 3.3b. We first convert 
the problem into an ordinary differential equation for si) Alex Ais 
is then solved iteratively using the method based on quasilinearization. 


The linearized problem at the n th iteration cycle is given by (3.79-81) 


with the following exceptions: 
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The superscript [n-1] in the last terms of (4.22) and (4.23) 
implies that the current iterate g in-1] is to be used for 6. 03 
and ¢, are constants defined in (3.73) and (3.75). The terms inside 
the curly brackets result from the unsteady term 96/at in equation (2.44). 


From the converged solution gti) we calculate the pressure distribution 


(FS) 


notisups sit mow -banissdo sd won nso. a ie 

tysvnoo t2vkt ow sdb lE% morse nt | oan ott 8 er 
otiT ox ink SPR eR fort sups (etn AeReeh venti Ws: + mide 
nottestsont fresvo no beasd beNtom ots onteu iapiastained: 
(f8-O1.£) yd nevip af efaya notdsistt Ay % ihe 35 metdong 


[r- gs AAT 


( &S.4) ; 2 


(€8.9)bns (SS.8)¥0 eirias desl “ond to 
ee 
53-6 tot beau. ad BY: ap blond 
_ a ae i 
ebtent amet siT .(2N.£) a CeA8 a 
7 yi 
(PBS) notteups nt t6\86 matey 


a. * 
notsudivgerb aweesxq oft ots Be. MEE rotay! 


95 


rae 

cl = j=2,9 (4.24) 
gti) 
J 

From this pti) we then calculate oe (j=2,J) using the finite dif- 


ference formulas given in Appendix C. 


4.2b Second Time Step 
At the second time step (i=2) the extrapolation formula 
(4.17) can not be used to approximate gl) | In order to maintain the 
same order of accuracy we use the following iterative scheme: 
p (2) and q (4) are obtained from (4.15) and (4.16) and a(2) is 


approximated by 


B81) (x) =p) (x) 
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The formula (4.25) has a truncation error of the order of 
— 2 2 
At, (a pyat’) (2) , eiigie! ) 


described in Section 4.2a is followed to obtain a better approximation 


With these values of a and ri the method 


for pt) and q(2), a (2) is then obtained from 
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which has the same order of truncation error as equation (4.17). The 
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calculations are repeated with these values of 
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The iteration is continued till pi?) converges within a 
predefined tolerance. Not more than three iterations was necessary to 


stabilize pi) to six significant digits. 


4.2c First Time Step 
At the first time step (i=1) we use the following iterative 
scheme: 
We start with the assumption pi) ¢x) = pee) (x), a") (x) = g6) (x), 
al) (,) = 0 and follow the method described in Section 4.2a to obtain 
a better approximation for pt!) and g(t), a1) (x) is then obtained 


from 


Ci) 7 x pl) X _-(0) x 
RN) = 48 ae eed (4.27) 


p 


which has a truncation error of order at, (apyat®) 1), We repeat the 
calculations with these values of mle gi!) and pl), 

The iteration is continued till pt!) converges within a 
predefined tolerance. To stabilize pl) to six significant digits, 


the number of iterations necessary was between four and six. 


4.3 Results of the Unsteady Problem 

The variation of p, f''(0), and H'(0) with time is presented 
in Figures 4.3 through 4.5. All curves are made for Pr = 0.72, y = 1.4, 
ue = Be eal and the change of angle from 2° to 6°. The initial 
and final states are also indicated. 


Figure 4.3 also indicates location of the condition t*ux/x* = ] 
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which corresponds to the discontinuities quoted in the literature [8,31]. 
The discontinuity is avoided by suitable formulation of the variables, 


and in particular by avoiding the variable of the type t*u*/x*. 
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CHAPTER V 
APPLICATION OF THE SOLUTIONS 


5.1 Aerodynamic Characteristics of Slender Wedge Wings 

Figure 5.1 shows the various forces acting on a slender wedge 
of semiwedge angle 8. at an angle of attack a. If the subscripts T 
and B represent the top and bottom surfaces of the wedge respectively, 


then the total normal force per unit width acting on the two sides are 


Figure 5.1 Aerodynamic Forces Acting on a Wedge 
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The total tangential force per unit width on the two sides are given 


where t* is the shear stress on the surface given by 


* 
q* = (u* ae (5.4) 
Mba 


Using the transformations of Chapter II we can write. 


Xi x Soha 
c* = 5 8 ie eG ; mp et | (5.5) 
a on. n=0 
Therefore from (5.3) we have 
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The lift and drag coefficients of the wedge wing can now be given as 
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Op 7 Oy + Oa (5.9a) 
6, = 8) - o (5.9b) 
Using (5.2) and (5.6) we have 
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where Ji> Jo represent the integrals 
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The lift to drag ratio is given by C /C)- 


Sor 


(T.2) 


(8.2) 


(s@.2) 


(d@.2) 


= aa t- 


™ "yas 


i > 

: * 
ee. 

i fs 7 
icf re 
. ec - = SON. 
oo + en 
- 


(dS? .2) 


104 
The position of the aerodynamic centre is given by 


L* L* 
nea b x* pk dx* - [> x* p& dx* 
eae os 
i 2 **00- 0 
or in the nondimensional form 
* = 
X = xac = J3\B Bir (5 14) 
a ik : Y Xpxl/2 : 
(J) pd 
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Here Js represents the integral 


(5.15) 


The rate of heat transfer to the wedge per unit area at any 
point of the surface is given by 


(5.16) 
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Making use of the transformations of Chapter II we can write 
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If Q* denote the total heat transfer rate per unit width to one side of 
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* — 3/2 
Diebae ede a eeu ear (5.18) 
F oxurL ZPY  &Q i 4°T ,B i 
where Jy represents the integral 
aH 
Felis ees ) dx (5.19) 
Jy | (x p on n=0 
fe) 
The integrals Jy> Jos Jz; Jy are of the form 
] 
W 
| x (x) dx (5.20) 
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Because of the special arrangement of the x-stations used in the nu- 
merical method, the function o(x) is given as a set of values - at 
the abscissae Xa j=1,2,.-.0. A modification of the Simpson's rule 


to evaluate the integrals of the form (5.20) is described in Appendix D. 
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APPENDIX A 
NUMERICAL METHODS FOR THE SERIES SOLUTION 


In this appendix we describe the numerical procedure for 


solving the coupled initial value problem of the form 
ts = F(f,F',f'' 5H.) (n)) (A.1a) 
TPR Sel Psi enlace ee) (A. 1b) 


with the initial values f(0), f'(0), f''(0), H(0O) and H'(0) specified. 
oy and bo represent known functions of n. A predictor-corrector method 
using the Falkner extrapolation formula for predictor and the Adams 
interpolation formula for corrector is employed. This scheme was 
originally chosen by Smith and Clutter [42] because of its suitability 
in high order problems of the form (A.1). 

We divide the interval n = 0 ton=n, by (N-1) uniform 


steps of size h. For any function $(n) we use the notation 
oy = o(n,) (A.2) 
where ip = (kel) Nerak= lees ose ; (A.3) 


Let us consider the general situation where the solution has 
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progressed up to Ne» k > 4. Extrapolating from the known values of f, 
H and their derivatives at the last four n-stations, we predict the un- 
known values at the station k+l. Representing the extrapolated values 


with the superscript ~ , the Falkner extrapolation formulas can be 


written as: 

Fedy = FED + gy (SSF 598) 43761 5-9F/ 4] (A.4a) 
Fi a \ h it tt 1 ! 
Fes = fr + aq LOSf, ~S9F) 0 peiten -9f L341 (a.4b) 

. ne 

Feed = fy + hf, + ==> hes [323F) ‘7 264F; | pr lsat, @ 27 38F | 31 (A.4c) 
Hea] = Hy + x [ 55H; K 99H,” 1*37Hy Oo 2-H) 13] (A.4d) 
Heep = Hie ay [55H'-59H! 1+37H! 5-9H! 4] (A.4e) 


The values of f''' and H'' at the station k+l can now be evaluated from 


A.1 and A.2 using the extrapolated values: 


Ue F(f k+]? flay feast k+] 01 Key) (A.5a) 
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The Adams interpolation formulas are now used to determine more accurate 


values of the unknowns at the station k+l: 
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fee 7 eh + gy DSF Oy! SELES (A.6a) 
Fee] ‘ bi cal oe CFP loFL | “Sip! it By ‘od (A. 6b) 
Fay = ft neL + day [38F LL FI7IFL 368 478 ,] (A.C) 
Hear = Hy + gy COAL, +19H! SHE +H] (A.6d) 
Hay = Hy + ag (Hey tIOHE-SHY HH 0] (A.6e) 


The values of Fay and He] are then calculated from A.1 and A.2. This 
completes the calculations at the station k+l and the solution pro- 
ceeds to the next station. 

The errors involved in the formulas (A.4) and (A.6) are of 
order (h°) or less. However, the errors in the extrapolation and the 
interpolation formulas are of opposite signs. Thus, the absolute 
value of the difference between the predicted and the corrected values 


provides an upper bound for the error at that particular step. If we 


define: 
E, = max {]$,-o, | sab = Tah ooh toate ees Hotas tine (A.7) 
then the error at any point in the solution is bounded by 


aay = max fe,» k=l, N} (A.8) 
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The choice of the step size h is governed by the maximum tolerance 
allowed one... In all the problems solved by this method, h = 2° 


(corresponding to N = 513 for tar 8) was found sufficient to bound 
Emax Within 1077, 

To start the solution by the predictor-corrector method 
described above, solutions at the stations k=1,2,3,4 are required. 
The values at station k=i (n=0) are provided by the initial conditions. 
The solutions at stations 2,3,4 are obtained by a Runge-Kutta method 
described by Romanelli [54]. Subroutine RKFA solves the initial value 


problem A.1 using the above algorithm. 
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APPENDIX B 
METHOD OF SOLVING THE LINEARIZED MOMENTUM 
AND ENERGY EQUATIONS 


Bl 


In this appendix we describe the finite difference method of 


solving the linearized problems (3.66,67) and (3.69,70). We use uni- 


form step size in the n-direction and follow the notation of equations 


(Ae 2,3)*. 
The following finite difference formulas are used: 


Five point symmetric formulas: 
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Five point unsymmetric formulas: 


(One forward and three backward points) 


: 3 2 
(LB) = (6, 4-66) p#12by 47100, 434,44)/2h> + 0(H*) 


UTE: 
(Lg) = (-0, thd, 9460, _1-200,4118,,,)/12h° + 0(h*) 
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Three point symmetric formulas: 
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Three point unsymmetric formulas: (Two forward points) 


), = (-39, +4645 ~b,49)/eh x 0(n°) 
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(B.2c) 
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The problem given by equations (3.66,67) can be written as: 


ae''' + b(nje'' + c(nje' + d(n)e = f(n) 
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e(0) = 0 (B.8a) 
e'(0) = 0 (B.8b) 
e'(n,) = 0 (B.8c) 


where a is a known constant and b,c,d,f are known functions of n. 


Replacing e''',e'',e' by the finite difference formulas we 
obtain the following system of linear algebraic equations for the vector 


e, (k=1 N) 
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(B.10) 


Zero elements of the matrix M are not shown. The elements in each 

row are normalized so that the magnitude of the largest element equals 1. 
The first, second and the last row of (B.9) are obtained 

from the boundary conditions (B.8) and the finite difference formulas 


(B.3) and (B.4). Thus we have 


A, ee B, = C, = D, a E = F = 0 (B.11a) 
= _ 49 ee ma ee Ae 2 
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The rows 3 to N-2 are obtained from the differential equation 


(oF) using the finite difference formulas (B.1). The elements are 


given by 
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D, = (-a + 5 hb, + £ hee.) /a (B. 12d) 

Eints (5 - by - hie) 0 (B.12e) 

F, = WF /a for k=3,4,...,Ne2 (B.12f) 
where a iS so chosen that 

max tA» Bie Cys Dy. E} = ] (B.13) 


For the (N-1)th row in (B.9) we have used the unsymmetric 


formulas (B.2). The elements are given by: 
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Where a is chosen to satisfy (B.13). 
The system of equations (B.9) is solved by Gaussian elimina- 


tion with partial pivoting. (B.10) is first transformed into 


Ue = F (B.15) 


where U is a NXN upper triangular matrix. Because of the pentadiagonal 
form of the matrix M, U is also of band structure with four superdiagonal 
elements in each row. From (B.15) e is calculated by back substitution. 
The procedure outlined above involving initial scaling of 
the elements and partial pivoting during elimination is recommended by 
Wilkinson [59], Forsythe and Moler [56] and is known to restrict the 
growth of round-off errors. The FORTRAN subroutine SLV5 was used to 
solve the linear system of the form (B.9). 
The linearized energy equation leads to the problem (3.69,70) 


which can be written as 
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g'(0) = 0 for insulated surface (B.17a) 

or g(0) = 0 for constant-temperature surface (Bel 7b) 
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Here a is a known constant and b,c,d are known functions of n. 
Since (B.16) is a second degree equation the three point 


formulas {3.5} are sufficient for truncation errors of 0(h*). The 
following system of linear algebraic equations is obtained for the 


vector ore Ck=1,N) 
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and Ay = OL By = ea Dread (B.19c) 


For k = 2,3,...,N-1 we have from (B.16) and (B.5) 


A= (a - Fb, )/a (B.20a) 

B, = (-2a + hc, )/a (B.20b) 

Cea Ob ia (B.20c) 

k Baek Tone 
D, = hod, fo (B. 20d) 
k k : 


The scaling factor a is so chosen that 
max {Ay ; By ; Cy 3 = ] (Bee) 


The system of equations (B.18) is solved by Gaussian elimina- 
tion with partial pivoting. The resulting upper triangular matrix 
is of band structure with two superdiagonal elements. The FORTRAN 
subroutine SLV3 was used to solve the system (B.18). Also the FORTRAN 
Subroutine MEEQN was used to solve the momentum and energy equations 


at any x-station using the method described in Section 3.3a. 
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APPENDIX C 
NUMERICAL METHODS FOR THE INTERACTION PROBLEM 


The problem given by equations (3.79,80) can be written as 


x 4 a(x)t = b(x) (C.1) 


t(x,) = 0 (C.2) 


where a(x) and b(x) are known functions of x. Using the arrangement 
of x stations described in (3.46), we solve the above problem by a 


finite difference method. The following formulas are used: 
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The coefficients Che (k,2 = 1,5) are functions of Ky and are 


given in Table C.1. The truncation errors of the formulas (C.3) are 
of 0(ax,"). 
Using the formulas (C.3) we obtain from (C.1) the following 


system of linear equations for T, estes. os 


ero iP T F, 
Bene ts ,| ES co 
fee, ce 0,| E, is F 
Ay Ba Cy Da Ey Ty Fy 
shy Ms ar ee aula ha 
Byes Pg ees eg dee ae, eal es 
Bjeoe eden Gasp Pgepeegeall 4| eo aie 
Aye Cj-eidel COMMIT | | eT Pet 
i (el pe Dae a iy F 

(C34) 


The elements of the matrix are given by 
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In equations (C.6-9) the quantities a are so chosen that 
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The matrix in the equation (C.4) is of the same form as that 
in equation (B.9). The linear system (B.9) is solved by the subroutine 
SLV5. The FORTRAN subroutine PRESS computes the pressure function p 
corresponding to a given I distribution by using the method described 


in Section 3.3b and this Appendix. 
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APPENDIX D 
EVALUATION OF THE INTEGRALS IN CHAPTER V 


In this appendix we develop a formula for the numerical 


evaluation of integrals of the type 


| x6(x) dx (D.1) 
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using a method similar to the Simpson's rule. The exponent w > - 1 


and os j=l,J are specified at a set of points Xs j=1,J as shown in 


Figure D.1. 
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Figure D.1 Set of Points Used to Evaluate Integral (D271) 
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For xi < xX eG we approximate ¢(x) by a second degree 


polynomial of the form 


o(x) = axe + bx +c (Daz) 


passing through the points (X5_ 905_4) (X35 505) and (X 54799541): 


Thus we have 
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The error involved in (0.3) is of the order ax? 9} | 
dx” j 
If J is odd, then the intervals between xy and x} are grouped 
two at a time and the formula (D.3) is applied on each group. Thus 


we have 
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Here the integral between x = 0 and x = xy has been calculated by 
using the same second degree curve that passes through the points 1,2,3. 
If J is odd, we use the above scheme for the integral up to 
X33 and calculate the remaining part by passing a third degree curve 
through the last four points. This gives the following formula: 
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The FORTRAN subroutine INTEG computes the integral (D.1) 


using this formula. 
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